# Set CRAN mirror to avoid errors
options(repos = c(CRAN = "https://cloud.r-project.org"))

# Function to safely install and load packages
install_and_load <- function(packages, bioconductor = FALSE) {
  for (package in packages) {
    # Check if package is installed
    if (!requireNamespace(package, quietly = TRUE)) {
      message(paste("Installing package:", package))
      
      # Install from appropriate source
      if (bioconductor) {
        if (!requireNamespace("BiocManager", quietly = TRUE)) {
          install.packages("BiocManager")
        }
        BiocManager::install(package, update = FALSE, ask = FALSE)
      } else {
        install.packages(package, repos = "https://cloud.r-project.org")
      }
    }
    
    # Load package
    suppressWarnings(suppressPackageStartupMessages(
      library(package, character.only = TRUE)
    ))
  }
}

# Install and load remotes package first
install_and_load("remotes")

# Regular CRAN packages
cran_packages <- c("survival", "lubridate", "ggfortify", "flexsurv", "remotes", 
                   "evaluate", "knitr", "dplyr", "ggplot2", "magrittr", "grid",
                   "broom.helpers", "data.table", "plotly", "DT", "flexdashboard",
                   "htmlwidgets", "visNetwork", "kableExtra", "RColorBrewer", "DT")

# Bioconductor packages
bioc_packages <- c("ggsurvfit", "gtsummary", "tidycmprsk", "ranger", 
                   "survminer", "survcomp", "maftools")

# Install and load packages
install_and_load(cran_packages)
install_and_load(bioc_packages, bioconductor = TRUE)

# Create bibliography directory if it doesn't exist
if (!dir.exists("bib")) {
  dir.create("bib")
}

# Generate package bibliography
all_packages <- c(.packages(), 'rmarkdown', 'knitr')
knitr::write_bib(x = all_packages, file = 'bib/packages.bib')

# Install maftools from GitHub if needed
if (!requireNamespace("maftools", quietly = TRUE)) {
  tryCatch({
    remotes::install_github("PoisonAlien/maftools", quiet = TRUE)
    library(maftools)
  }, error = function(e) {
    message("Couldn't install maftools from GitHub, trying Bioconductor version instead")
  })
}

# Define file paths
maf_file_path <- "/Users/sultanalharbi/Library/CloudStorage/OneDrive-Personal/Projects/Thesis_Chapters/Chapter 1 (cfDNA-TACE for HCC)/TACE_Tables/HCC_TACE_maftools.maf"
clinical_data_path <- "/Users/sultanalharbi/Library/CloudStorage/OneDrive-Personal/Projects/Thesis_Chapters/Chapter 1 (cfDNA-TACE for HCC)/TACE_Tables/ClinicData_Maf_Final.tsv"

# Check if files exist
maf_exists <- file.exists(maf_file_path)
clinical_exists <- file.exists(clinical_data_path)

if (!maf_exists || !clinical_exists) {
  warning("One or more files not found. Please check the paths.")
  if (!maf_exists) warning(paste("MAF file not found:", maf_file_path))
  if (!clinical_exists) warning(paste("Clinical data file not found:", clinical_data_path))
} else {
  # Read MAF file with all mutation types
  HCC_Final <- read.maf(
    maf = maf_file_path,
    clinicalData = clinical_data_path,
    rmFlags = FALSE,
    useAll = TRUE,
    verbose = TRUE,
    isTCGA = FALSE,
    vc_nonSyn = c(
      "Missense_Mutation", "5'UTR", "3'UTR", "5'Flank", "In_Frame_Del",
      "Frame_Shift_Del", "Frame_Shift_Ins", "Splice_Site", "Nonsense_Mutation", 
      "Intron", "Silent"
    )
  )
  
  # Read MAF file without intron/silent mutations
  HCC_Final_non_intron <- read.maf(
    maf = maf_file_path,
    clinicalData = clinical_data_path,
    rmFlags = FALSE,
    useAll = TRUE,
    verbose = TRUE,
    isTCGA = FALSE,
    vc_nonSyn = c(
      "Missense_Mutation", "5'UTR", "3'UTR", "5'Flank", "In_Frame_Del",
      "Frame_Shift_Del", "Frame_Shift_Ins", "Splice_Site", "Nonsense_Mutation"
    )
  )
  
  # Read MAF file with default non-synonymous definitions
  HCC_Final_non_synonymous <- read.maf(
    maf = maf_file_path,
    clinicalData = clinical_data_path,
    rmFlags = FALSE,
    useAll = TRUE,
    verbose = TRUE,
    isTCGA = FALSE,
    vc_nonSyn = NULL
  )
  
  # Load clinical data
  ClinicalData <- read.delim(clinical_data_path)
  ClinicalData <- as.data.table(ClinicalData)
  
  # Associate clinical data with MAF object
  HCC_Final@clinical.data <- ClinicalData
  HCC_Final_non_intron@clinical.data <- ClinicalData
  HCC_Final_non_synonymous@clinical.data <- ClinicalData
  
  # Print summary of loaded data
  cat("Data successfully imported:\n")
  cat(paste0("- Total samples: ", length(unique(HCC_Final@clinical.data$Tumor_Sample_Barcode)), "\n"))
  cat(paste0("- Total patients: ", length(unique(HCC_Final@clinical.data$Pateint)), "\n"))
  cat(paste0("- Total mutations (all types): ", nrow(HCC_Final@data), "\n"))
  cat(paste0("- Total non-intron mutations: ", nrow(HCC_Final_non_intron@data), "\n"))
  cat(paste0("- Total non-synonymous mutations: ", nrow(HCC_Final_non_synonymous@data), "\n"))
  
  # List of HCC genes of interest
  HCC_genes <- c(
    "APC", "ARID1A", "ARID2", "ATM", "AXIN1", "CDKN2A", "CTNNB1", "HNF1A", 
    "JAK1", "KEAP1", "LZTR1", "PIK3CA", "PTEN", "RB1", "SF3B1", 
    "SMARCA4", "TERT", "TP53"
  )
if (!require("remotes")) install.packages("remotes")
remotes::install_github("uham-bio/UHHformats", build_vignettes = TRUE)
  # Define color palette for mutation types
  fabcolors <- RColorBrewer::brewer.pal(n = 12, name = 'Set3')
  names(fabcolors) <- c(
    "Missense_Mutation", "5'UTR", "3'UTR", "5'Flank", "In_Frame_Del",
    "Frame_Shift_Del", "Frame_Shift_Ins", "Splice_Site", "Intron", 
    "Nonsense_Mutation", "Silent"
  )
  fabcolors <- list(Variant_Classification = fabcolors)
}

#==============================================================================
# FUNCTION DEFINITIONS
#==============================================================================
#' Create MAF subsets based on clinical queries
#' 
#' @param maf MAF object to subset
#' @param subset_definitions Named list of query expressions
#' @return None (assigns objects to global environment)
create_named_subsets <- function(maf, subset_definitions) {
  for (name in names(subset_definitions)) {
    assign(name, subsetMaf(maf = maf, clinQuery = subset_definitions[[name]]), envir = .GlobalEnv)
  }
}

#' Calculate comprehensive statistics for MAF subsets
#' 
#' @param maf_objects List of MAF objects
#' @param reference_maf Reference MAF for percentage calculations
#' @param category_map Mapping of subset names to categories
#' @return Data frame with statistics for each subset
calculate_maf_stats <- function(maf_objects, reference_maf = NULL, category_map = NULL) {
  # Initialize an empty list to store individual data frames
  stats_list <- list()
  
  # Get reference statistics if available
  ref_samples <- NA
  ref_mutations <- NA
  if(!is.null(reference_maf) && inherits(reference_maf, "MAF")) {
    ref_samples <- length(unique(reference_maf@clinical.data$Tumor_Sample_Barcode))
    ref_mutations <- nrow(reference_maf@data)
  }
  
  # Process each subset
  for(name in names(maf_objects)) {
    tryCatch({
      maf_obj <- maf_objects[[name]]
      
      # Get category if provided
      category <- "Other"
      if(!is.null(category_map) && name %in% names(category_map)) {
        category <- category_map[[name]]
      }
      
      # Calculate basic statistics
      samples <- unique(maf_obj@clinical.data$Tumor_Sample_Barcode)
      n_samples <- length(samples)
      n_mutations <- nrow(maf_obj@data)
      
      # Skip if no data
      if(n_samples == 0 || n_mutations == 0) {
        stats_list[[name]] <- tibble(
          Category = category,
          Subset = name,
          Samples = n_samples,
          Mutations = n_mutations,
          MedianMutsPerSample = 0,
          TopGene = "None",
          TopGeneMutationCount = 0,
          SamplePct = ifelse(!is.na(ref_samples), round(n_samples/ref_samples*100, 1), NA),
          MutationPct = ifelse(!is.na(ref_mutations), round(n_mutations/ref_mutations*100, 1), NA)
        )
        next
      }
      
      # Calculate mutation metrics
      muts_per_sample <- table(maf_obj@data$Tumor_Sample_Barcode)
      median_muts <- median(as.numeric(muts_per_sample))
      
      # Find top gene
      gene_counts <- sort(table(maf_obj@data$Hugo_Symbol), decreasing = TRUE)
      top_gene <- names(gene_counts)[1]
      top_gene_count <- gene_counts[1]
      
      # Calculate percentages
      sample_pct <- ifelse(!is.na(ref_samples), round(n_samples/ref_samples*100, 1), NA)
      mut_pct <- ifelse(!is.na(ref_mutations), round(n_mutations/ref_mutations*100, 1), NA)
      
      # Create tibble for this subset (tibbles don't have row names)
      stats_list[[name]] <- tibble(
        Category = category,
        Subset = name,
        Samples = n_samples,
        Mutations = n_mutations,
        MedianMutsPerSample = median_muts,
        TopGene = top_gene,
        TopGeneMutationCount = top_gene_count,
        SamplePct = sample_pct,
        MutationPct = mut_pct
      )
    }, error = function(e) {
      warning(paste("Error processing", name, ":", e$message))
    })
  }
  
  # Combine using bind_rows from dplyr to avoid row names issues
  if(length(stats_list) > 0) {
    stats_df <- bind_rows(stats_list)
  } else {
    stats_df <- tibble(
      Category = character(),
      Subset = character(),
      Samples = integer(),
      Mutations = integer(),
      MedianMutsPerSample = numeric(),
      TopGene = character(),
      TopGeneMutationCount = integer(),
      SamplePct = numeric(),
      MutationPct = numeric()
    )
  }
  
  return(stats_df)
}

#==============================================================================
# SUBSET DEFINITIONS
#==============================================================================
# Define cohort categories
subset_categories <- list(
  "Treatment" = c("pre_tace_samples", "post_tace_samples", "progression_samples"),
  "Etiology" = c("MASLD_samples", "HCV_samples", "HBV_samples", "Alcohol_samples", 
                 "ALLsample_Viralsamples", "ALLsample_NoViral", "Alcohol_sample_MASLD_samples",
                 "ALLsample_No_MASLD", "ALLsample_No_Alcohol", "ALLsample_No_HBC"),
  "BCLC" = c("BCLC_A", "BCLC_B", "BCLC_C", "ALLBCLC_NoC", "ALLBCLC_NoA"),
  "Survival" = c("Patient_OS_survival", "Patient_PFS_survival")
)

# Define subset query expressions
subset_definitions <- list(
  # Treatment condition subsets
  pre_tace_samples = "Condition1 == 'Pre_TACE'",
  post_tace_samples = "Condition1 == 'Post_TACE'",
  progression_samples = "Condition1 == 'progression'",
  
  # Disease etiology subsets
  MASLD_samples = "Aetiology == 'MASLD'",
  HCV_samples = "Aetiology == 'HCV'",
  HBV_samples = "Aetiology == 'HBV'",
  Alcohol_samples = "Aetiology == 'Alcohol'",
  ALLsample_Viralsamples = "Aetiology %in% c('HCV','HBV')",
  ALLsample_NoViral = "Aetiology %in% c('Alcohol','MASLD','Other')",
  Alcohol_sample_MASLD_samples = "Aetiology %in% c('Alcohol','MASLD')",
  ALLsample_No_MASLD = "Aetiology %in% c('Alcohol','HCV','HBV','Other')",
  ALLsample_No_Alcohol = "Aetiology %in% c('MASLD','HCV','HBV','Other')",
  ALLsample_No_HBC = "Aetiology %in% c('MASLD','HCV','Alcohol','Other')",
  
  # BCLC stage subsets
  BCLC_A = "BCLC %in% c('A')",
  BCLC_B = "BCLC %in% c('B')",
  BCLC_C = "BCLC %in% c('C')",
  ALLBCLC_NoC = "BCLC %in% c('A','B')",
  ALLBCLC_NoA = "BCLC %in% c('B','C')",
  
  # Survival status subsets
  Patient_OS_survival = "OS_Status == '1'",
  Patient_PFS_survival = "PFS_Status == '1'"
)

#==============================================================================
# CREATE COHORT SUBSETS
#==============================================================================
# Create all subsets at once
cat("### Creating MAF subsets\n")

create_named_subsets(HCC_Final, subset_definitions)

cat("✓ Created", length(subset_definitions), "MAF subsets\n\n")

# Create category mapping
category_map <- list()
for(cat in names(subset_categories)) {
  for(subset in subset_categories[[cat]]) {
    category_map[[subset]] <- cat
  }
}

# Get all MAF objects
maf_objects <- list()
for(subset_name in unlist(subset_categories)) {
  if(exists(subset_name) && inherits(get(subset_name), "MAF")) {
    maf_objects[[subset_name]] <- get(subset_name)
  }
}

#==============================================================================
# CALCULATE STATISTICS
#==============================================================================
cat("### Calculating statistics\n")

subset_stats <- calculate_maf_stats(maf_objects, reference_maf = HCC_Final, category_map = category_map)
# Ensure it's a tibble to avoid row name issues
subset_stats <- as_tibble(subset_stats)
cat("✓ Calculated statistics for", nrow(subset_stats), "subsets\n\n")

#==============================================================================
# GENERATE SUMMARY REPORT
#==============================================================================
# Overall dataset summary
cat("# Hepatocellular Carcinoma Cohort Statistics\n\n")

# Reference dataset information
ref_samples <- length(unique(HCC_Final@clinical.data$Tumor_Sample_Barcode))
ref_mutations <- nrow(HCC_Final@data)
ref_median_muts <- median(table(HCC_Final@data$Tumor_Sample_Barcode))

cat("<div style='background-color: #f8f9fa; padding: 15px; border-radius: 5px; margin-bottom: 20px;'>")

cat("<strong>Reference Dataset Summary:</strong><br>")

cat(ref_samples, "samples with", ref_mutations, "mutations<br>")

cat("Median mutations per sample:", ref_median_muts, "<br>")

cat("</div>\n\n")

#==============================================================================
# CATEGORY-SPECIFIC STATISTICS
#==============================================================================
# Create tables and visualizations by category
for(cat in unique(subset_stats$Category)) {
  cat(paste0("## ", cat, " Cohorts\n\n"))
  
  # Filter statistics for this category
  cat_stats_filtered <- subset_stats %>% 
    filter(Category == cat) %>%
    arrange(desc(Samples))
  
  # Create a completely fresh tibble without any possibility of row names
  cat_stats <- tibble(
    Subset = cat_stats_filtered$Subset,
    Samples = cat_stats_filtered$Samples,
    Mutations = cat_stats_filtered$Mutations,
    MedianMutsPerSample = cat_stats_filtered$MedianMutsPerSample,
    TopGene = cat_stats_filtered$TopGene,
    TopGeneMutationCount = cat_stats_filtered$TopGeneMutationCount,
    SamplePct = cat_stats_filtered$SamplePct,
    MutationPct = cat_stats_filtered$MutationPct
  )
  
  cat("### Summary Table\n")
  # Add row.names=FALSE parameter here
  knitr::kable(cat_stats, 
               caption = paste(cat, "Cohort Statistics"),
               format = "html",
               row.names = FALSE) %>%  # Explicitly set row.names=FALSE
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                 full_width = FALSE) %>%
    column_spec(1, bold = TRUE) %>%
    print()
  
  cat("\n### Visualizations\n")
  # Sample count visualization
  p1 <- ggplot(cat_stats, aes(x = reorder(Subset, -Samples), y = Samples, fill = Subset)) +
    geom_bar(stat = "identity") +
    geom_text(aes(label = Samples), vjust = -0.5, size = 3.5) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = "none") +
    labs(title = paste("Sample Counts for", cat, "Cohorts"),
         x = NULL, y = "Number of Samples") +
    scale_fill_brewer(palette = "Set2")
  
  # Mutation rate visualization  
  p2 <- ggplot(cat_stats, aes(x = reorder(Subset, -MedianMutsPerSample), 
                            y = MedianMutsPerSample, fill = Subset)) +
    geom_bar(stat = "identity") +
    geom_text(aes(label = MedianMutsPerSample), vjust = -0.5, size = 3.5) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1),
          legend.position = "none") +
    labs(title = paste("Mutation Rate for", cat, "Cohorts"),
         x = NULL, y = "Median Mutations per Sample") +
    scale_fill_brewer(palette = "Set2")
  
  # Display plots
  print(p1)
  print(p2)
  
  cat("\n\n")
}

#==============================================================================
# CROSS-CATEGORY COMPARISONS
#==============================================================================
cat("## Cross-Category Comparisons\n\n")

# Sample size comparison across categories
# Use proper tibble-compatible approach for factor reordering
subset_stats_ordered <- subset_stats %>%
  mutate(Subset = factor(Subset, levels = Subset[order(Samples, decreasing = TRUE)]))

# Overall sample size comparison
cat("### Overall Sample Count Comparison\n")

overall_plot <- ggplot(subset_stats_ordered, aes(x = Subset, y = Samples, fill = Category)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Samples), vjust = -0.3, size = 3) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Cohort Sizes Across All Categories",
       x = NULL, y = "Number of Samples") +
  scale_fill_brewer(palette = "Set2")

print(overall_plot)

# Mutation rate comparison
cat("\n### Overall Mutation Rate Comparison\n")

mutation_plot <- ggplot(subset_stats, aes(x = reorder(Subset, -MedianMutsPerSample), 
                                        y = MedianMutsPerSample, fill = Category)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = MedianMutsPerSample), vjust = -0.3, size = 3) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Mutation Rates Across All Categories",
       x = NULL, y = "Median Mutations per Sample") +
  scale_fill_brewer(palette = "Set2")

print(mutation_plot)

#==============================================================================
# INTERACTIVE DATA EXPLORATION
#==============================================================================
cat("\n## Interactive Data Table\n\n")

if(requireNamespace("DT", quietly = TRUE)) {
  # Create a completely new tibble from scratch to avoid inheriting row names
  display_stats <- tibble(
    Category = subset_stats$Category,
    Subset = subset_stats$Subset,
    Samples = subset_stats$Samples,
    Mutations = subset_stats$Mutations,
    MedianMutsPerSample = subset_stats$MedianMutsPerSample,
    TopGene = subset_stats$TopGene,
    TopGeneMutationCount = subset_stats$TopGeneMutationCount,
    SamplePct = subset_stats$SamplePct,
    MutationPct = subset_stats$MutationPct
  )
  
  DT::datatable(display_stats,
                caption = "Comprehensive Cohort Statistics",
                filter = "top",
                options = list(
                  pageLength = 20,
                  dom = 'Blfrtip',
                  buttons = c('copy', 'csv'),
                  lengthMenu = list(c(10, 20, -1), c('10', '20', 'All'))
                ),
                rownames = FALSE) %>%
    DT::formatStyle(
      'Samples',
      background = DT::styleColorBar(c(0, max(display_stats$Samples)), 'steelblue'),
      backgroundSize = '98% 88%',
      backgroundRepeat = 'no-repeat',
      backgroundPosition = 'center'
    ) %>%
    DT::formatStyle(
      'Mutations',
      background = DT::styleColorBar(c(0, max(display_stats$Mutations)), 'forestgreen'),
      backgroundSize = '98% 88%',
      backgroundRepeat = 'no-repeat',
      backgroundPosition = 'center'
    ) %>%
    DT::formatRound(columns = c('MedianMutsPerSample', 'SamplePct', 'MutationPct'), digits = 1)
} else {
  cat("*Note: Install the DT package for interactive data tables*")
}

# Create a dashboard with key metrics
valueBox_html <- function(value, caption, icon, color) {
  paste0(
    '<div class="section-box" style="background-color:', 
    switch(color,
           "primary" = "#e8f0f9",
           "info" = "#e3f2fd",
           "success" = "#e8f8ef",
           "warning" = "#fff8e8",
           "#e8f0f9"), 
    '; padding: 15px; border-radius: 5px; margin: 10px; flex: 1; min-width: 250px; box-shadow: 0 1px 3px rgba(0,0,0,0.12);">',
    '<div style="display:flex; align-items:center;">',
    '<div style="font-size:2.5em; margin-right:15px; color:', 
    switch(color,
           "primary" = "#3498db",
           "info" = "#0D47A1",
           "success" = "#27ae60",
           "warning" = "#f39c12",
           "#3498db"), 
    ';">',
    '<i class="fa ', icon, '"></i>',
    '</div>',
    '<div>',
    '<div style="font-size:1.8em; font-weight:bold;">', value, '</div>',
    '<div style="font-size:1.2em; color:#7f8c8d;">', caption, '</div>',
    '</div>',
    '</div>',
    '</div>'
  )
}

if (exists("HCC_Final") && exists("HCC_Final_non_intron") && exists("HCC_Final_non_synonymous")) {
  # Number of samples
  num_samples <- length(unique(HCC_Final@clinical.data$Tumor_Sample_Barcode))
  # Number of mutations
  num_mutations <- nrow(HCC_Final@data)
  # Median mutations per sample
  median_muts <- median(HCC_Final@variants.per.sample$Variants)
  # Number of samples non-intron
  num_samplesnon_intron <- length(unique(HCC_Final_non_intron@clinical.data$Tumor_Sample_Barcode))
  # Number of mutations non-intron
  num_mutationsnon_intron <- nrow(HCC_Final_non_intron@data)
  # Median mutations per sample non-intron
  median_mutsnon_intron <- median(HCC_Final_non_intron@variants.per.sample$Variants)
  # Number of samples non-synonymous
  num_samplesnon_synonymous <- length(unique(HCC_Final_non_synonymous@clinical.data$Tumor_Sample_Barcode))
  # Number of mutations non-synonymous
  num_mutationsnon_synonymous <- nrow(HCC_Final_non_synonymous@data)
  # Median mutations per sample non-synonymous
  median_mutsnon_synonymous <- median(HCC_Final_non_synonymous@variants.per.sample$Variants)
} else {
  cat("Data not properly loaded. Please check the data import step.")
}

if (requireNamespace("gtsummary", quietly = TRUE) && requireNamespace("DT", quietly = TRUE) && exists("HCC_Final")) {
  library(gtsummary)
  library(dplyr)
  library(DT)
  
  # Get all unique patients from the first Pre_TACE sample for each patient
  unique_patients <- HCC_Final@clinical.data %>%
    as.data.frame() %>%
    filter(Condition1 == "Pre_TACE") %>%
    distinct(Pateint, .keep_all = TRUE)
  
  # If not all patients have Pre_TACE samples, add their first available sample
  if (nrow(unique_patients) < length(unique(HCC_Final@clinical.data$Pateint))) {
    missing_patients <- setdiff(unique(HCC_Final@clinical.data$Pateint), unique_patients$Pateint)
    
    for (patient in missing_patients) {
      patient_samples <- HCC_Final@clinical.data[HCC_Final@clinical.data$Pateint == patient, ]
      unique_patients <- rbind(unique_patients, patient_samples[1, ])
    }
  }
  
  # Select relevant columns for summary
  clin_cols <- c("Age", "Gender", "Aetiology", "Viral_Not_Viral", "Cirrhosis", 
                 "CTP", "BCLC", "Tumour_Size", "Response", "OS_time",
                 "OS_Status", "PFS_time", "PFS_Status", "PS", "HAP", "PVT")
  
  # Create summary table for available columns
  available_cols <- intersect(clin_cols, colnames(unique_patients))
  
  if (length(available_cols) > 0) {
    # Create categorization for variables
    var_labels <- list(
      Age = "Age (years)",
      Gender = "Sex",
      Aetiology = "Etiology",
      Viral_Not_Viral = "Viral Status",
      Cirrhosis = "Cirrhosis",
      CTP = "Child-Turcotte-Pugh Class",
      BCLC = "Barcelona Clinic Liver Cancer Stage",
      Tumour_Size = "Tumor Size (cm)",
      Response = "Treatment Response",
      OS_time = "Overall Survival (months)",
      OS_Status = "Overall Survival Status",
      PFS_time = "Progression-Free Survival (months)",
      PFS_Status = "Progression Status",
      PS = "Performance Status",
      HAP = "Hepatoma Arterial-embolization Prognostic Score",
      PVT = "Portal Vein Thrombosis"
    )
    
    # Subset to available variables
    var_labels <- var_labels[names(var_labels) %in% available_cols]
    
    # Create summary statistics in an organized way
    tbl_summary_data <- unique_patients %>%
      select(all_of(available_cols)) %>%
      tbl_summary(
        by = NULL,
        label = var_labels,
        statistic = list(
          all_continuous() ~ "{median} ({p25}, {p75})",
          all_categorical() ~ "{n} ({p}%)"
        ),
        digits = all_continuous() ~ 1,
        missing = "no"
      ) %>%
      modify_header(label = "**Characteristic**", stat_0 = "**Value (N={N})**") %>%
      bold_labels() %>%
      as_gt() %>%
      gt::tab_header(title = "Clinical Characteristics of HCC Patients Treated with TACE")
    
    # Print the table
    tbl_summary_data
    
    # Create interactive table
    DT::datatable(
      unique_patients %>% select(all_of(available_cols)),
      options = list(
        pageLength = 10,
        scrollX = TRUE,
        dom = 'Blfrtip',
        buttons = c('copy', 'csv', 'excel'),
        columnDefs = list(
          list(className = 'dt-center', targets = "_all")
        )
      ),
      filter = 'top',
      class = 'cell-border stripe hover',
      caption = htmltools::tags$caption(
        style = 'caption-side: top; text-align: center; font-size: 16px; font-weight: bold; margin-top: 20px;',
        'Interactive Table: HCC Patient Data'
      )
    )
  } else {
    cat("Clinical data columns not found in the expected format.")
  }
} else {
  cat("Required packages (gtsummary, DT) not available or data not properly loaded.")
}

if (exists("HCC_Final_non_intron") && exists("HCC_genes")) {
  # Generate oncoplot for HCC genes
  tryCatch({
    oncoplot(
      maf = HCC_Final_non_intron, 
      genes = HCC_genes,
      clinicalFeatures = c('Aetiology', 'BCLC', 'Cirrhosis'),
      sortByAnnotation = TRUE,
      annotationColor = NULL,
      removeNonMutated = FALSE,
      showTumorSampleBarcodes = TRUE,
      fontSize = 0.7,
      SampleNamefontSize = 0.8,
      legendFontSize = 1,
      annotationFontSize = 0.8,
      showTitle = TRUE,
      titleText = "Mutation Profile of Key HCC Genes",
      showPct = TRUE
    )
  }, error = function(e) {
    cat("Error in generating oncoplot:", e$message, "\n")
    cat("Trying with fewer options.\n")
    
    # Simplified oncoplot
    oncoplot(
      maf = HCC_Final_non_intron,
      genes = HCC_genes,
      showTumorSampleBarcodes = TRUE
    )
  })
}

# ------------------------------------------------------------------
# Extended oncoplot – top 18 genes with detailed clinical tracks
# ------------------------------------------------------------------
if (exists("HCC_Final_non_intron")) {

  #--- make sure the clinical annotation names really exist -----------
  anno_cols <- colnames(HCC_Final_non_intron@clinical.data)

  # Use the long name only if it is available, otherwise fall back to BCLC
  clin_feats <- c(
    'Aetiology',
    'Condition1',
    if ('Barcelona_Clinic_Liver_Cancer' %in% anno_cols) {
      'Barcelona_Clinic_Liver_Cancer'
    } else if ('BCLC' %in% anno_cols) {
      'BCLC'
    } else {
      NULL
    }
  )

  #--------------------------------------------------------------------
  tryCatch({

    oncoplot(
      maf                 = HCC_Final_non_intron,
      top                 = 18,
      genes               = c("APC", "ARID1A", "ARID2", "ATM", "AXIN1", "CDKN2A",
                              "CTNNB1", "HNF1A", "JAK1", "KEAP1", "LZTR1",
                              "PIK3CA", "PTEN", "RB1", "SF3B1", "SMARCA4",
                              "TERT", "TP53"),
      altered             = FALSE,
      drawRowBar          = TRUE,
      drawColBar          = TRUE,
      includeColBarCN     = TRUE,
      clinicalFeatures    = clin_feats,
      annotationColor     = NULL,
      draw_titv           = FALSE,
      showTumorSampleBarcodes = TRUE,
      barcodeSrt          = 90,
      gene_mar            = 5,
      anno_height         = 1,
      legend_height       = 4,
      sortByAnnotation    = FALSE,
      groupAnnotationBySize = TRUE,
      sortByMutation      = FALSE,
      keepGeneOrder       = FALSE,
      GeneOrderSort       = TRUE,
      removeNonMutated    = TRUE,
      bgCol               = "#CCCCCC",
      borderCol           = "white",
      fontSize            = 0.7,
      SampleNamefontSize  = 0.9,
      titleFontSize       = 1.5,
      legendFontSize      = 1.2,
      annotationFontSize  = 1.2,
      showTitle           = TRUE,
      titleText           = "Top 18 Mutated Genes – Detailed Oncoplot",
      showPct             = TRUE
    )

  }, error = function(e) {
    message("Oncoplot (top 18) could not be generated: ", e$message)
  })
}

laml.titv = titv(maf = HCC_Final_non_intron, plot = TRUE, useSyn = TRUE)

#plot titv summary
plotTiTv(
  res = laml.titv,
  plotType = "both",
  sampleOrder = HCC_Final_non_intron@clinical.data$Tumor_Sample_Barcode,
  color = NULL,
  showBarcodes = TRUE,
  textSize = 0.8,
  baseFontSize = 1,
  axisTextSize = c(1, 1),
  plotNotch = FALSE
)

png("/Users/sultanalharbi/Library/CloudStorage/OneDrive-Personal/Projects/Thesis_Chapters/Chapter 1 (cfDNA-TACE for HCC)/TACE_Results/somatic_interactions.png", width = 12, height = 9, units = "in", res = 300)
par(mar = c(5, 8, 4, 2))  # Increase left margin (2nd value)

somaticInteractions(
  maf = HCC_Final_non_intron,
  genes = HCC_genes,
  pvalue = c(0.05, 0.1),
  leftMar = 8,
  topMar = 8,# Default is 4, increase to 8 or more
  fontSize = 1,
  showSigSymbols = TRUE,
  returnAll = FALSE  # Make sure this is FALSE to get a plot instead of data
)

dev.off()

# Define the list of HCC genes to include
HCC_genes <- c(
  "APC", "ARID1A", "ARID2", "ATM", "AXIN1", "CDKN2A", "CTNNB1", "HNF1A", 
  "JAK1", "KEAP1", "LZTR1", "PIK3CA", "PTEN", "RB1", "SF3B1", 
  "SMARCA4", "TERT", "TP53"
)

# Define the extract_vaf_data function
extract_vaf_data <- function(maf_object, group_label) {
  # Extract the variants data from the MAF object
  variants <- maf_object@data
  
  # Select relevant columns and add group label
  result <- variants %>%
    select(Hugo_Symbol, tumor_f) %>%
    mutate(Group = group_label)
  
  return(result)
}

# Extract and process data as before
pre_data <- extract_vaf_data(pre_tace_samples, "Pre-TACE")
post_data <- extract_vaf_data(post_tace_samples, "Post-TACE")

# Combine the data
combined_data <- bind_rows(pre_data, post_data)

# Filter out NA values and restrict to only the specified HCC genes
combined_data <- combined_data %>%
  filter(!is.na(tumor_f)) %>%
  filter(Hugo_Symbol != "NA" & !is.na(Hugo_Symbol) & Hugo_Symbol != "") %>%
  filter(Hugo_Symbol %in% HCC_genes)  # Add this line to filter for HCC genes only

# Important: Set factor levels to control the order
combined_data$Group <- factor(combined_data$Group, levels = c("Pre-TACE", "Post-TACE"))
# Set gene order to match the original vector order
combined_data$Hugo_Symbol <- factor(combined_data$Hugo_Symbol, levels = HCC_genes)

# Create side-by-side boxplots with Pre-TACE first
vaf_plot_ordered <- ggplot(combined_data, aes(x = Hugo_Symbol, y = tumor_f, fill = Group)) +
  geom_boxplot(position = position_dodge(width = 0.8), width = 0.7, outlier.size = 0.5, na.rm = TRUE) +
  scale_fill_manual(values = c("Pre-TACE" = "#3498db", "Post-TACE" = "#e74c3c")) +
  labs(
    title = "VAF Comparison in HCC-Related Genes",
    x = "Gene",
    y = "Variant Allele Frequencies (VAF)",
    fill = ""
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
    axis.text.y = element_text(size = 10),
    legend.position = "top",
    legend.title = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# Save the plot
ggsave("/Users/sultanalharbi/Library/CloudStorage/OneDrive-Personal/Projects/Thesis_Chapters/Chapter 1 (cfDNA-TACE for HCC)/TACE_Results/pre_post_vaf_comparison_hcc_genes.png", vaf_plot_ordered, width = 10, height = 6)

# For the stacked vertical panels version, the same factor setting ensures correct order
vaf_plot_stacked <- ggplot(combined_data, aes(x = Hugo_Symbol, y = tumor_f, fill = Group)) +
  geom_boxplot(outlier.size = 0.5) +
  scale_fill_manual(values = c("Pre-TACE" = "#3498db", "Post-TACE" = "#e74c3c")) +
  facet_grid(Group ~ ., scales = "free_y") +  # Pre-TACE will now appear on top
  labs(
    title = "VAF Comparison in HCC-Related Genes",
    x = "Gene",
    y = "Variant Allele Frequencies (VAF)",
    fill = ""
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
    axis.text.y = element_text(size = 10),
    strip.text = element_text(size = 12, face = "bold"),
    legend.position = "none",
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

# Save the plots
ggsave("/Users/sultanalharbi/Library/CloudStorage/OneDrive-Personal/Projects/Thesis_Chapters/Chapter 1 (cfDNA-TACE for HCC)/TACE_Results/vaf_comparison_stacked_panels_hcc_genes.png", vaf_plot_stacked, width = 10, height = 8)
print(vaf_plot_ordered)